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Abstract 

A viscous- in vise id interactive coupling method is 
used for the computation of unsteady transonic flows in- 
volving separation and reattachment. A lag-entrainment 
integral boundary layer method is used with the tran- 
sonic small disturbance potential equation in the CAP- 
TSDV code. Efficient and robust computations of steady 
and unsteady separated flows, including steady separa- 
tion bubbles and self-excited shock-induced oscillations 
are presented. The buffet onset boundary for the NACA 
0012 airfoil is accurately predicted and shown computa- 
tionally to be a Hopf bifurcation. Shock-induced oscilla- 
tions are also presented for the 1 8 percent circular arc air- 
foil. The oscillation onset boundaries and frequencies are 
accurately predicted, as is the experimentally observed 
hysteresis of the oscillations with Mach number. This 
latter stability boundary is identified as a jump phenom- 
enon. Transonic wing flutter boundaries are also shown 
for a thin swept wing and for a typical business jet wing, 
illustrating viscous effects on flutter and the effect of sep- 
aration onset on the wing response at flutter. Calculations 
for both wings show limit cycle oscillations at transonic 
speeds in the vicinity of minimum flutter speed indices. 

Introduction 

The onset of transonic shock-induced flow separa- 
tion is known to be associated with a variety of nonclas- 
sical aeroelastic instability and response phenomena, 1-3 
referred to variously as; single degree of freedom flutter, 
limited-amplitude flutter, limit cycle oscillations (LCO), 
control surface buzz, and buffeting (onset). A charac- 
teristic of the "instabilities" involved is a tendency to 
grow to a constant or bounded "limit amplitude" which 
can vary from a nuisance level to levels large enough 
to cause structural failure. In the latter case, the non- 
classical response, genetically referred to herein as LCO, 
is typically observed near the flutter boundary, making 
a distinction between the two response mechanisms dif- 
ficult. Edwards 4,5 reviewed these features of transonic 
aeroelasticity, concluding that i.) computational capabil- 
ity for such cases would require modeling of dynamically 
separating and reattaching viscous boundary layers and 
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ii.) such capability was not yet mature for wings or more 
complete configurations. 

Transonic separation can also lead to the occurrence 
of self-excited shock-induced oscillations (SIO) for per- 
fectly rigid structures; well documented in the case of 
two-dimensional (2-D) airfoils. 6 ' 10 The unsteady airloads 
generated are quite large, occur for narrow ranges of 
transonic Mach number (and angle of attack), and have 
characteristic frequencies which can be near those in- 
volved in flutter. The ability to calculate SIO condi- 
tions may be necessary in order to treat the LCO phe- 
nomena. Refs. 10-12 report calculations of SIO using 
2-D unsteady Navier-Stokes codes. While shock oscilla- 
tions were obtained, the calculations were sensitive to de- 
tailed modeling issues (wind tunnel geometry, grid spac- 
ing, etc.), computed frequencies tended to be low, and 
results were given only for isolated conditions within the 
interior of SIO regions. More recent calculations of SIO 
for these two airfoils are given in Refs. 13-15 where 
very good agreement with SIO onset boundaries and fre- 
quencies are obtained. 



Fig. 1 Sketch of shock-boundary layer interaction. 

Interactive Boundary Layer Modeling (IBLM) pro- 
vides an alternative to such direct computation of flows 
involving viscous shear layers. Separate computations 
are made for an inner viscous boundary layer region and 
an outer inviscid flow region as illustrated in Fig. 1 . Sub- 
script “e” denotes the “edge” of the boundary layer, while 
superscripts “i” and “v” denote inviscid and viscous vari- 
ables. Ref. 16 developed an integral boundary layer lag- 
entrainment method to compute displacement thickness 
<5* which was used to update the flow tangency boundary 
condition of the inviscid solver. This "direct" solution 
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method for the entrainment equation becomes singular at 
flow separation and "inverse" computation methods 17 " 22 
have been developed in attempts to treat flow separa- 
tion. Carter 17 introduced an iterative, relaxation method 
for coupling the inner and outer region edge velocities 
for steady flows. Melnik and Brook 18 incorporated the 
free shear layer closure modeling of Ref. 19 into their 
steady inverse lag-entrainment method, using Carter’s 17 
coupling method to interact with a steady full potential 
code. Results 18 indicate that only small regions of sep- 
arated flow could be treated. LeBalleur and Girodroux- 
Lavigne 19 developed a semi-implicit method for coupling 
unsteady boundary layer and Transonic Small Distur- 
bance (TSD) potential equations, showing examples of 
transonic separated flow including SIO conditions. Un- 
steady lag-entrainment boundary layer and TSD equations 
were also coupled in Refs. 20 and 21 using a "quasi- 
simultaneous" method. Although cases involving sepa- 
rated flow on oscillating airfoils are reported, no exam- 
ples of SIO calculations are given. Both of these un- 
steady interactive methods 19, 21 rely upon explicit bound- 
ary condition treatment available during the final z -sweep 
of the 2-D Alternating-Direction Implicit (ADI) solution 
algorithm. Locally linear relations between the inner 
and outer flow variables are developed and pointwise 
iterations with relaxation are used. Bartels 14, 15 has re- 
cently developed an IBLM with a fully unsteady finite- 
difference boundary layer model and presents many SIO 
calculations. 

References 23-27 present further examples of steady 
IBLM coupled with steady potential and Euler equation 
solvers. The 2-D transonic separated flow cases shown 
involve separation onset; shock-induced separation bub- 
bles and trailing-edge separation. Comparisons 23 24 in- 
dicate the IBLM solutions are quite competitive with 
Navier-Stokes solutions. 


cous flow, whose coupling requires active control ele- 
ments in order to minimize the coupling error between 
the two systems. Reference 32 gives details of this de- 
velopment and calculations of two types of self-excited 
shock oscillations about rigid airfoils. This paper summa- 
rizes these SIO calculations for the NACA 0012 and for 
the 18 percent thick circular arc airfoils. The extention 
of the method implemented in the three-dimensional code 
in a stripwise fashion and termed the CAP-TSDV code, 
is described below. Calculations of wing flutter are pre- 
sented for two cases: a 4 percent thick wing flutter model, 
and a thicker typical business jet wing flutter model. 


Mathematical Method 

Details of the inviscid flow equation, the boundary 
layer equations, the modifications to boundary conditions, 
and the IBLM coupling procedure are summarized in this 
section. Further details are given in Ref. 32. 


Inviscid Flow Equation 

The CAP-TSD potential equation code is used in 
this analysis. The CAP-TSD code uses an approximate 
factorization algorithm for time-accurate solution of the 
unsteady TSD equation 29 


dfo , + d Jl JL d Jl = o 

dt dx dy dz 


( 1 ) 


where <f> is the inviscid-disturbance velocity potential 
and 


/o = -A<f>t - B<j> x 

( 2 a) 

fi = E\<j> x + Fi4>l + 

( 2 b) 

fl = <t>y E\<j>x<t>y 

( 2 c) 

h = 4>z 

( 2 d) 


Howlett 28 developed an IBLM coupled with the 
CAP-TSD 29 (Computational Aeroelasticity Program - 
Transonic Small Disturbance) potential equation code. 
Howlett’s modification of Rizzetta’s 30 quasi-steady, di- 
rect, lag-entrainment IBLM was incorporated in the CAP- 
TSD code in a stripwise fashion. Further modifications 31 
involved an inverse solution procedure using Carter’s 
coupling method. Experience in calculating unsteady sep- 
arated transonic flows was similar to that noted above for 
steady conditions; only small amounts of flow separation 
could be treated without code failure. 

Edwards 32 reported a new coupling method, using 
Howlett's inverse solution procedure, capable of treat- 
ing transonic SIO conditions for airfoils. The new cou- 
pling method was based upon the observation that at the 
transonic flow conditions of interest, the flowfield is fre- 
quently inherently unsteady with oscillating shocks and 
dynamically separating and reattaching boundary layers. 
The IBLM is thus regarded as a simulation of two dy- 
namic systems, the outer inviscid flow and the inner vis- 


where A = M 2 , B = 2 A/ 2 , E\ = 1 — M 2 , F\ — 
— i[3-(2-7 )M 2 ]M 2 , (?! = — 5 M 2 , and H 1 = 
— M 2 . For the 2-D version of the code, fa = G\ = 
0. The code contains modifications to these coefficients 
developed by Batina 33 to approximate the effects of shock 
generated entropy and vorticity. 


The boundary conditions on the airfoil/wing 
wake are 

and 

tf=St+Sf ; 

x le < X < X te> 2 = 0* 

(3) 

A^>* — 0 \ 

X > X te , Z = 0* 

(4) 

A {4>x + <f>t) = 0 ; 

X > X te . X — 0* 

(5) 

where the superscript ± refers to the airfoil upper 

and 


lower surface, S(x,t) ( S(x,y,t )) denotes the airfoil 
(wing) shape, and A(- • •) indicates the jump in (• • •) 
across the wake. Nonreflecting far-field boundary con- 
ditions are also used. There have been extensive applica- 
tions of the 3-D code for unsteady aerodynamic and flutter 
analysis of wings and complete configurations. 29, 35 
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Lag-Entrainment Boundary Layer Equations 

The effect of a turbulent viscous boundary layer is 
modeled in the quasi-steady manner of Green et al. 16 by 
solving a set of ordinary differential equations in x for the 
integral boundary layer quantities: momentum thickness 
0; shape factor 77; and entrainment coefficient Ce 




dx 


2\/U* 

xx, t 


( 6 ) 


= ( c * - \ H ' C ) + 


,dCE 

dx 


= f {hTh;1 {Ct)eqo 


(±*£\ 

dx ) EQ 


,2 1 + 


- a cy z ]+ 


. 

2 


i + 


( 8 ) 


Equation (8) is taken from Ref. 19 and differs slightly 
from that given in Ref. 16. The subscript e refers 
to quantities at the edge of the boundary layer and the 
inviscid surface velocity gradient 4>xx is obtained from 
the outer flow solver. The various closure parameters in 
these equations are given in Ref. 31. Equations (6-8) are 
suitable for attached flow boundary layers and provide the 
boundary layer displacement thickness 


6* = HO 


(9) 


This provides a "direct" calculation of the viscous mod- 
ification to the airfoil shape to be implemented in the 
boundary conditions discussed below. 


Viscous Boundary Conditions and Wake 

The coupling between the inner and outer solutions 
is through the boundary conditions on the airfoil/wing and 
wake. Equations (3) and (4) are modified as follows: 

<j>f = St + St + s; ; Xie < X < Xu, z = 0 ± (10) 
A<t>, = A{6 m x ) ; x>x u , 2 = 0 ± (11) 


Viscous Equations for Separated Flow 

Correlations of the IBLM equations given above 
with experiment indicate the onset of flow separation oc- 
curs at conditions where H ~ 2.0 — 2.3 and C j ~^0 . At 
these values of 77, Eqs. (6-8) become singular (dH/dH i 
becomes infinite), and alternative solution procedures are 
necessary. The present study utilizes the "inverse IBLM 
as implemented by Howlett, 31 which closely follows that 
of Melnik and Brook. 18 Eqs. (6-7) are inverted to provide 
solutions for 77* and u v tt as functions of the boundary 


layer displacement thickness, represented by a perturba- 
tion mass flow parameter fn = p* t xi\6*. Equation (8) 
remains unchanged, giving symbolically 


du\ „ _ 1 dm 

-& = Fl + F2 7fi~fc 
dH_ 1 dm 

— _ F 3 -f d x 

^ = I(f 5 + f 6 *U) 


( 12 ) 

(13) 

(14) 


The Fi variables arc nonlinear functions of the param- 
eters modeling the closure conditions for attached and 
separated flow and are defined in Ref. 31. For separated 
flows, these closure conditions are based upon the work 
of Melnik and Brook, 18 which closely follows the de- 
velopment of LeBalleur and Girodroux-Lavigne, 19 where 
additional details may be found. 


Numerical Implementation 

From the leading edge of the airfoil, the boundary 
layer is approximated by the turbulent boundary layer on 
a flat plate. At a specified point, numerical integration of 
Eqs. (12-14) is implemented with a fourth-order Runge- 
Kutta method. For the Mach number range studied, it 
was found that the inverse equations, in conjunction with 
the coupling method described below, converged rapidly 
for attached flow upstream of regions of flow separation 
(and also for downstream regions of reattached flow). 
This obviated the use of Eqs. (6-8) thus circumventing 
the numerically troublesome switching between the direct 
and inverse equations in separating flow regions, where 
the largest parameter gradients occur. 

Interactive Boundary Layer Coupling Method 

Since the intended applications of the IBLM include 
cases of wing flutter, including SIO and LCO, the cou- 
pling method was developed based on the observation that 
at the transonic flow conditions of interest, the flowfield 
is frequently inherently unsteady, displaying oscillating 
shocks and separating and reattaching boundary layers. 
The interacting boundary layer method is thus regarded 
as a simulation of two dynamic systems, the outer in- 
viscid flow and the inner viscous flow, whose coupling 
requires active control elements in order to minimize the 
coupling error between the two systems. The elements 
utilized, illustrated in analog block diagram fashion in 
Fig. 2, include a variable gain integral control element 
for the displacement thickness and a first order smoothing 
filter for the momentum thickness estimate. The coupling 
equations are thus 


e{x) - u*(i) - u\(x) = «"(*) - ( 15 ) 

= K(Ki(6')e (16) 

at 
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Equations (17-19) are treated as ordinary differential 
equations in time at each span wise location and are im- 
plemented using digital filtering methods. The nonlinear 
gain, K X (S *), is scheduled on the local displacement 
thickness in order to enhance the abilitiy of the coupling 
method to follow moving “wedge-shaped” separated flow 
regions. Smooth blending of the attached and separated 
flow closure relations over the range 1.5 < H <2.5 is 
implemented via linear interpolation. Finally, the dis- 
placement thickness in the wake was modeled to have 
exponential decay. Numerical experiments 32 indicated 
relative insensitivity of SIO solutions to the parameters 
of the wake modeling. 



Fig. 2 Schematic diagram of variable gain, integral control, 
viscous-invisdd interative coupling method. 

For the 3-D CAP-TSDV code, Eqs. 12-14 and 
15-17 are solved independently at each spanwise chord 
station on the wing. This is accomplished at each time 
step, within the Newton linearization iteration loop of the 
approximate factorization solution algorithm of the CAP- 
TSD code 36 . 

Model Descriptions and Results 
Airfoil Models 

Calculations with the variable gain integral control 
coupling method are shown for two airfoils and two wing 
flutter models. The NACA 0012 airfoil and the 1 8% thick 
circular arc airfoil demonstrate cases of steady shock- 
induced separation bubbles, steady trailing-edge separa- 
tion, and SIO for rigid airfoils. Extensive experimental 
results have been published for both airfoils and particular 
attention is drawn to the reports of McDevitt and Okuno 8 
and McDevitt. 6,7 For both airfoils, the wind tunnel walls 
were contoured in order to approximate steady free-air 
streamline flow (wall boundary layer mass removal using 
suction was also utilized for the NACA 0012 tests) and 
regions of Mach number and angle of attack for the onset 
and quenching of SIO are reported. 


The airfoil calculations were obtained on a 212 x 100 
grid in x— z space. The grid extends ±20 chords in x and 
±25 chords in z. On the airfoils, 146 grid points were 
distributed in order to capture the large shock excursions 
involved in SIO: aft of midchord a uniform spacing of 
0.5 percent chord was used; forward of 40 percent chord 
the uniform spacing was l| percent chord. 

NACA 0012 Airfoil Calculations 

Reference 8 reports results from tests of the NACA 
0012 airfoil for Mach numbers from approximately 0.7 
to 0.8 and at angles of up to six degrees, sufficient to 
induce buffet onset. Reynolds number based on chord, 
Re Cf ranged from 1 to 14 million. The test section uti- 
lized flexible upper and lower walls, including boundary 
layer suction, in order to minimize wall interference and 
provide contouring compatible with free air streamlines. 
The onset of unsteady "buffet" flow is shown to occur 
along a well defined boundary of Mach number vs. angle 
of attack. The pressure waveforms in the buffet region 
were erratic at low Reynolds number, but for the higher 
Reynolds numbers well developed cyclic shock oscilla- 
tions were seen. 

Calculations have been made with the 2-D version 
of CAP-TSDV for steady conditions just below the buffet 
onset boundary and for SIO conditions at and well beyond 
the boundary. Results were obtained for time steps of 
At = 0.05 and for Re c = 10 million. Unless otherwise 
noted, the inverse IBLM was initiated at 10 percent chord, 
Kt = 0.00010, and 5 Newton iterations were performed. 
Figure 3a compares calculated and experimental surface 
pressures for M = 0.775 and a = 2.05°. (Data Set 5 
of Ref. 8). The agreement is very good, capturing the 
pressure levels, shock location, and separation bubble 
at the foot of the shock. Calculated boundary layer 
parameters, Fig. 3b, indicate the rapid growth of 6* 
through the shock and the effect of the separation bubble 
is clearly seen in C/. 

Calculations shown in Fig. 4 indicate a well-defined 
buffet onset boundary in excellent agreement with the 
experimental boundary. 8 The present results are an im- 
provement over those of Ref. 37, which indicated onset 
to occur approximately 0.5 degrees higher than experi- 
ment. Also, the CAP — TSDV code is able to compute 
the onset boundary for 0°, where the shock excursion 
amplitude is largest. 

Figure 5 shows the development of the SIO from 
“zero” flow field conditions for M = 0.775 and a = 
4.0°. While a dominant SIO frequency of k b = 0.25 
(based on semichord) can be seen, the lift coefficient 
exhibits significant nonlinearity. 

The development of the SIO with increasing angle 
of attack shows a Hopf bifurcation onset characteristic, 
illustrated in Fig. 6. The peak-to-peak oscillation ampli- 
tude of the lift coefficient bifurcates at a = 2.3°. The 
amplitude grows quickly to a maximum of 0.3 while the 
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SIO frequency remains rather constant Further penetra- 
tion of the buffet region, to a = 4.0 deg., leads to in- 
creasing frequency, in good agreement with experiment 8 . 
The bifurcation amplitude curve is reversible, with the 
same amplitudes being obtained for decreasing angles. 


method must thus be evaluated based on agreement with 
experimental results and other calculations (see Ref. 32 
for more details). Figures 4 and 6 indicate that the present 
method gives good agreement with experimental buffet 
onset and frequency. This is very encouaging, since it 
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Fig. 3 Comparison of calculated and experimental results 

for NACA 0012 airfoil at M = 0.775, a = 2.05°, Rec = 10 7 . 

For conditions below the buffet onset boundary, 
fully converged viscous solutions are obtained, includ- 
ing stable shock-induced separation bubbles as in Fig. 
3. Since the boundary layer equations are quasi-steady, 
while the CAP — TSDV code is time accurate, it is not 
surprising that fully converged solutions cannot be ob- 
tained at each time step for conditions above the buffet 
onset boundary. The capability of the current coupling 


Fig. 4 Comparison of calculated and experimental buffet 
onset boundaries for the NACA 0012 airfoil, 
is just such difficult computational conditions for which 
the coupling method is intended. Note particularly the 
ability of the method to achieve a significant penetration 
into the buffet region as seen in Fig. 6. 
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Fig. 5 Lift and moment coefficient time histories for the 
NACA 0012 airfoil at SIO conditions: M = 0.775, or = 4.0 , 
Rec = 10 7 . 
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Fig. 6 Reduced frequency and amplitude of lift coefficient 
versus angle. NACA 0012 airfoil, M = 0.775, Re« = 10 7 . 

Eighteen Percent Thick Circular 
Arc Airfoil Calculations 

Numerous wind tunnel tests have reported SIO for 
the 18 percent thick circular arc airfoil. References 6 
and 7 document SIO for 0.73 < M < 0.78 and for 
Reynolds numbers from 1 to 17 million. Tests were con- 
ducted in a blow-down wind tunnel with walls contoured 
to simulate free air streamline conditions. The oscilla- 
tions were at a characteristic frequency of kb % 0.47, 6 
varying little with Mach number, and decreasing slightly 
for nonzero angles. A sensitivity of the range of SIO 
was found for Re c « 1 — 6 million. A significant fea- 
ture is a hysteresis of the SIO phenomenon with regard 
to increasing or decreasing Mach numbers. A number 
of SIO calculations for this airfoil are discussed by Ed- 
wards and Thomas. 11 Successful SIO calculations for iso- 
lated Mach numbers have been made using Navier-Stokes 
codes 11, 12,38 and a TSD code with an IBLM. 15 These 
Navier-Stokes results were all low in calculated SIO fre- 
quency, kt « 0.40, while the TSD result was k h « 0.34. 
More recent calculations 13, 14 show very good agreement 
with experiment. 

Calculations have been made with the CAP-TSDV 
2-D code for conditions encompassing the experimental 
SIO region up to M = 0.78. Effects of Mach number, 
integral coupling gain, grid size, wake modeling, and 
disturbance amplitude upon the SIO have been studied 32 . 


Unless otherwise noted, dt = 0.025, Re c = 10 million, 
K 6 = 0.00015, the IBLM was initiated at 10 percent 
chord, and 5 Newton iterations were performed at each 
time step. 




Fig. 7 Lift and moment coefficient time histories for the 
18% thick circular arc airfoil at SIO conditions: M = 0.76, 
Re c = 10 7 . 



Fig. 8 Comparison of calculated and experimental SIO 
reduced frequencies for the 18% thick circular arc airfoil. 

The buildup of the SIO at M = 0.76 from ''zero” 
flow field conditions is shown in Fig. 7. The buildup of 
the SIO to its limit-amplitude occurs in three shock oscil- 
lation cycles over 2.5 chordlengths of travel. Comparison 
of Figures 4 and 7 indicates more nearly sinusoidal air- 


6 


loads for this case, involving upper and lower surface 
shocks oscillating antisymmetrically, than for the NACA 
0012 airfoil case, which involves a single upper surface 
shock motion. 

Calculated and experimental 6 reduced frequencies 
for these SIO are compared in Fig. 8. The calculated 
frequencies are in very good agreement with experiment, 
including the mild decrease in frequency with increasing 
Mach number. 

The amplitude of the shock excursion is approxi- 
mately 40 percent chord with the aftmost shock posi- 
tion increasing about 5 percent chord from M = 0.74 
to 0.78. The displacement thickness at the trailing edge 
varies from 1.5—11 percent chord. All of the computed 
shocks are "strong" shocks as discussed by McDevitt 6 
(supersonic-to-subsonic flow) with nearly constant sub- 
sonic pressure levels aft of the shocks to the trailing- 
edge. The forward motion of the shock for all three 
Mach numbers stalls just forward of the crest of the airfoil 
at midchord, coinciding with rising pressures aft of the 
shock and formation of a suction pressure region, lead- 
ing to formation of the new shock. The computed shock 
speeds during their forward movement are quite constant, 
with speed increasing slightly with increasing Mach num- 
ber. The increasing shock speed with Mach number, in 
conjunction with the larger shock excursion mentioned 
above leads to almost constant periods for the SIO. The 
gain used for these results, K& = 0.00015, gives very 
good agreement 32 with the experimental SIO frequency 
of Jfc = 0.47. 

Wind tunnel tests 6 ’ 7 show an interesting hysteresis 
of the onset of the SIO phenomenon depending upon in- 
creasing or decreasing Mach number. The new coupling 
method is able to reproduce this effect. Figure 9 super- 
imposes calculated SIO onset and quenching boundaries 
with those from Ref. 6. For increasing Mach number at 
Re c = 10 million, SIO onset occurs for M = 0.76 exper- 
imentally and for M = 0.755 computationally. For de- 
creasing Mach number, the SIO quenches for M ~ 0.73 
experimentally and for M = 0.735 computationally. For 
lower Reynolds numbers, experimental results indicate a 
narrowing of the SIO region, presumably due to transi- 
tion effects. The present calculations show an indication 
of such an effect of Reynolds number, although smaller 
in magnitude. No attempt has been made to modify the 
present method with transition modeling features. 

Comparison of the SIO onset boundaries in Figs. 
4 and 9 indicates that different mechanisms are operat- 
ing in the two cases. Inspection of steady flow condi- 
tions occurring for increasing Mach numbers within the 
hysteresis region, 0.73 < M < 0.76, confirms this to 
be the case. Figure 10 gives steady computed surface 
pressure and boundary layer parameters for M = 0.74. 
Trailing-edge separation at x/c = 0.85 is clearly evi- 
dent. Also evident is the near separation at x/c = 0.60 
induced by the shock-boundary layer interaction. This 


pattern of incipient shock-induced and trailing-edge flow 
separation offers an explanation of the SIO hysteresis 
for this airfoil. With increasing Mach number, separa- 
tion initiates at the trailing-edge while the weak shock 
on the airfoil is not yet strong enough to separate the 
boundary layer. Investigation of boundary layer parame- 
ters (not shown) indicates the local minimum in the skin 
friction coefficient, C/, observed just downstream of the 
shock decreases continuously with increasing Mach num- 
ber. When this local minimum C/ « 0 at M = 0.755, the 
spontaneous SIO discussed above develop. With the SIO 
established, oscillations persist for decreasing Mach num- 
ber until M = 0.735 where the minimum value (steady) 
of Cj at the shock is « 0.0010. 


no SIO SIO 


O ♦ present calc. 



(a) increasing Mach number. 


no SIO SIO 

O ♦ present calc. 



(b) decreasing Mach number. 

Fig. 9 Regions of SIO for increasing and decreasing Mach 
number for the 18% thick circular arc airfoil. 
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The hysteresis of the SIO for 0.73 < M < 0.76, 
coupled with the flow separation patterns discussed above 
suggests that oscillations might be induced, within this re- 
gion, by appropriate perturbations. Figure 1 1 shows that 
this is indeed the case. For M = 0.74, perturbations 
were introduced by simulating trailing-edge control sur- 
face motions, 6j (exponentially shaped pulses were used 
for the l/4c flap). For 6j < 10°, stable decaying airload 
oscillations are seen, whereas the oscillations quickly lock 
onto the SIO waveform for 6f > 10°. Thus an amplitude 
threshold, jump phenomenon is identified as one of the 
SIO mechanisms for this airfoil. 



(a) pressure coefficient. 




(b) boundary layer variables. 

Fig. 10 Steady calculated results for the 18% thick circular 
arc airfoil within SIO hysteresis region, M = 0.74, Rec = 10 7 . 


d • dog. 




Fig. 11 Lift and moment time histories for the 18% thick 
circular arc airfoil exhibiting jump phenomenon due to flap 
pulses, M = 0.74. 


Wing Flutter Models 

The first wing flutter model, shown in Fig. 12, is the 
AGARD Standard Aeroelastic Configuration 39 * 40 which 
was tested in the Transonic Dynamics Tunnel (TDT) at 
NASA Langley Research Center. It is a semispan wall- 
mounted model having a quarter-chord sweep angle of 
45 deg., a panel aspect ratio of 1.65, and a taper ratio of 
0.66. The wing had a NACA 65A004 airfoil sction and 
was constructed of laminated mahogany. The wing is 
modeled structurally using the first four natural vibration 
modes, with natural frequencies ranging from 9.6 Hz for 
the first bending mode to 91.54 Hz for the second torsion 
mode The CAP-TSDV calculations were performed on a 
150 x 30 x 80 point computational grid with 100 points 
along each of 15 spanwise chords on the wing. Other 
computational conditions were: nondimensional time step 
dt = 0.05, one Newton iteration, and K b = 0.00030. 

The second wing flutter model, shown in Fig. 13, is 
a typical business jet configuration also tested in the TDT. 
The semispan wing-fuselage model was mounted on the 
wind tunnel sidewall and tested in air, with experimental 
flutter data obtained for Mach numbers from 0.628 to 
0.888. The wing has a taper ratio of 0.29 and a midchord 
sweep of 23 degrees. The airfoil thickness varies from 
13 percent at the symmetry plane (for the extended wing- 
alone configuration analyzed) to 8.5 percent at the wing 
tip. Six natural vibration modes were included in the 
calculations, with frequencies ranging from 4.3 Hz to 62.7 
Hz. The CAP-TSDV calculations were performed on a 
100 x 50 x 80 point computational grid with 45 points 
along each of 33 spanwise chords on the wing. Other 
computational conditions were: nondimensional time step 
dt = 0.03, one Newton iteration, and K h = 0.00010. 
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Fig. 12 Plan view of AGARD Wing 445.6 Standard Aeroe- 
lastic Configuration. 



Fig. 13 Business jet flutter model mounted in NASA 
Langley Transonic Dynamics Tunnel. 


AGARD Wing 445.6 Flutter Calculations 

Model tested in air. The majority of published 
calculations for this model (actually a series of models 
with similar planforms) are for the “weakened model #3“ 
tested in air, since this test covered the largest transonic 
speed range and showed a significant transonic dip ef- 
fect. Figure 14 gives new results from CAP-TSDV (large 
square symbols) and includes results from other NASA 
Langley Research Center studies 41,42 for comparison. It 
is informative to discuss these results with respect to two 
Mach number ranges. 

— M < 1.0. The CAP-TSDV results for these four 
Mach numbers are in excellent agreement with experi- 
ment. The comparisons with other calculations illustrate 
details which appear to be relevant to such aeroelastic 
analysis (below Mach one) in general. As a point of ref- 


Q_ experiment 
CAP-TSD, linear 41 



(a) flutter speed index. 



Mach Number 

(b) frequency. 

Fig. 14 Comparison between experimental and calculated 
flutter speed index and frequency for the AGARD Wing 
445.6 tested in air. 

erence, linear theory flutter calculations for this case 41 
are also in very good ageement with the experimental 
points except for that at M — 0.96 where the calcu- 
lation is unconservative. This is to be expected, since 
nonlinear transonic effects for this thin wing are confined 
to Mach numbers very near one. For the three lower 
Mach numbers, the CAP-TSDV and CFLSD-Euler 42 re- 
sults are in excellent agreement. The inviscid calcula- 
tions for M = 0.901 and 0.96 increasingly overpredict 
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the drop in the transonc flutter boundary. At M = 0.96 
the two viscous calculations, from CAP-TSDV and the 
CFL3D-Navier-Stokes, indicate that viscous modeling is 
required to correct this overprediction. This effect ap- 
pears to grow in importance for thicker wings, as will be 
seen below. 

— M > 1.0. Fewer calculations have been pub- 
lished dealing with this very low supersonic Mach num- 
ber range and the results summarized in figure 14 in- 
dicate that more work is needed to understand differ- 
ences between experiment and calculation. Again, linear 
CAP-TSD results are in very good agreement with experi- 
ment. Higher level calculations have been less successful 
for these two low supersonic cases. The inviscid Euler 
equation results 42 are very unconservative and the vis- 
cous Navier-Stokes result corrects only 80 percent of the 
discrepancy at M = 1.141 42 . Invisid CAP-TSD results 
(not shown) are somewhat higher in flutter speed index 
than the Euler results and the effect of including the vis- 
cous boundary layer at M = 1.141 is less effective in 
correcting the discrepancy than that shown in Figure 14. 
A number of factors may be considered in discussing this 
discrepancy. The flutter boundary for this model is quite 
sensitive to Mach number here, as noted by the steep gra- 
dient seen in these two flutter points. In addition, none 
of the codes attempt to model the details of the cut-off 
tip of the model. Finally, for these very low supersonic 
Mach numbers, wind tunnel interference effects are not 
well understood. 

Model tested in heavy gas. Since the above results 
for the model tested in air resulted in somewhat unrealis- 
tically large mass ratios and small reduced flutter frequen- 
cies, it was desirable to obtain results for the weakened 
models #5 and #6” which were tested in heavy gas and 
had more reasonable ranges of mass ratio and frequency. 
CAP-TSDV calculations for these cases are shown in 
Fig. 15. Again, for these cases with M < 1.0, the 
CAP-TSDV results are in excellent agreement with ex- 
periment for M = 0.74 and 0.92. Due to issues discussed 
above for very low supersonic Mach numbers, calculation 
have not been attempted for the third experimental Mach 
number of 1.0. Instead calculations at M = 0.94, 0.95, 
and 0.96 revealed an interesting minimum feature in the 
flutter speed index parameter at M = 0.95. Further nu- 
merical experimentation at M = 0.96 revealed nonlinear 
response features. It was found that the estimated damp- 
ing of the flutter mode was dependent upon amplitude. 
Figure 16 shows a simulated wing “tip rap” response for 
M = 0.96 and Q = 0.75 psi. The interpolated flut- 
ter dynamic pressure from the experimental data for this 
Mach number is Q * = 0.757 psi. The early portion 
of the response indicates positive damping of the flut- 
ter mode and a higher frequency mode. The damping 
of the flutter mode decreases as the response amplitude 
decays to approximately 0.12 inches peak-to-peak, where 
stable limit cycle oscillations persist. This limit cycle be- 
havior was further studied by sequentially increasing the 


dynamic pressure between computed runs from Q — 0.5 
to 0.81 psi.. The resulting tip deflection time history is 



Mach number 


(a) flutter speed index. 



(b) frequency. 

Fig. 15 Comparison between experimental and calculated 
flutter speed index and frequency for the AGARD Wing 
445.6 tested in heavy gas. 

shown in Figure 16. Eleven computer runs with a total of 
22,000 time steps were calculated. The dynamic pressure 
was incremented as indicated in steps between restarted 
runs. For Q < 0.60 psi. the response is damped and 
for Q = 0.70 psi. small neutrally stable oscillations are 
seen. With Q increased to 0.78 psi. slowly divergent 
oscillations develop and with further increase to 0.81 psi. 
the divergent oscillations grow with increased negative 
damping until the amplitude reaches approximately 0.12 
inches peak-to-peak. The growth of the oscillations then 
quenches and it appears that a limit cycle condition will 
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again develop, although further calculations are needed 
to fully establish this feature. 



Fig. 16 Calculated AGARD Wing 445.6 tip response in heavy gas for M = 0.96 and Q = 0.75 psi. 



Fig. 17 Calculated AGARD Wing 445.6 response in heavy gas for M = 0.96 and increasing dynamic pressure. 
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This limit cycle behavior for this model was only 
observed for the highest calculated Mach number, M = 

0.96 which lies on the “backside” of the small transonic 
dip seen in Fig. 15. At this Mach number and for the 
wing motions calculated, the flow is fully attached with 
no significant transonic features. The boundary layer cou- 
pling method performed well, with well-converged dis- 
placement thickness profiles. Numerical flow visualiza- 
tions of the wing pressure surfaced details which are pos- 
sibly key to this nonlinear response behavior. At this 
Mach number and for this thin wing significant regions 
of near sonic flow develop adjacent to the wing upper 
and lower surfaces as the wing oscillates. Very high fre- 
quency upstream moving pressure waves are seen in the 
visualizations which are consistent with forward propa- 
gating Mach waves. At a given point on the wing, the 
frequency of these pressure waves is 10-20 times the 
flutter mode frequency for this case. The amplitudes of 
these calculated limit cycles is small and and no mention 
of such behavior is reported 40 . It is unlikely that such 
small motions, even if present, would have been detected 
since they would have been heavily masked by the model 
response to tunnel turbulence. 

Business Jet Wing Flutter Calculations 

The business jet wing flutter model shown in Fig- 
ure 4 was tested in the Transonic Dynamics Tunnel at 
NASA Langley Research Center. Gibbons 43 presents flut- 
ter calculations for the model including spatial and tempo- 
ral convergence studies, and surface pressure coefficient 
comparisons for rigid and statically deformed cases, using 
TSD, Euler, and Navier-Stokes methods. For the present 
study, the effect of including viscous effects using the 
CAP-TSDV code was investigated. 

The model was constructed from aluminum plate 
with fiberglass wrapped foam providing the airfoil con- 
tour. The wing was mounted low on the side-wall 
mounted fuselage model which had a circular cross- 
section with a conical aft end. The wing root angle- 
of-attack was varied during the test to minimize loading. 
The maximum angle needed for this purpose was 0.2 de- 
grees at the highest tested Mach number. This root angle 
was used for the calculations described below. This re- 
sulted in calculated static tip deflections of -1.33 in. at 
M = 0.628 and -hi. 35 in. at M — 0.888. The Reynolds 
numbers for these two Mach numbers were 2.17 million 
and 1.14 million respectively, based on the 2.0 ft. root 
chord. The model had a 4.4 ft. semispan. 


Contour plots of the upper and lower wing surface 
pressure, displacement thickness, and skin friction are 
shown in Fig. 18 for M = 0.888. Note the lower surface 
leading edge suction peak and mild inboard shock seen in 
the Fig. 18a. The lower surface displacement thickness 
is similar to the upper surface with maximum thicknesses 
below one percent except near the root where the lower 
surface shock produces a thickness of approximately 1.5 
percent root chord. The skin friction in Fig. 18c. reflects 
these features seen in the displacement thickness and is 
informative regarding closeness to separation. The lower 
surface trailing-edge is separated at the root and there 
is small separation bubble just inboard of the tip and 
aft of the leading-edge suction peak. The skin friction 
coefficient is low in the trailing-edge region of the upper 
surface, reaching a minimum near 88 percent span. This 
region and the lower surface separation bubble are key 
in the effect of amplitude upon flutter mode response 
described below. 

Calculated flutter speed indices and frequencies ver- 
sus Mach number are compared with experiment in Fig. 
19. The linear CAP-TSD, Euler, and Navier-Stokes re- 
sultes are from Gibbons 43 while the four CAP-TSDV data 
points are new. Comparison of these flutter boundaries 
leads to similar observations as for the 445.6 wing: 

1. Inviscid calculations agree among themselves 
and are in very good agreement with experiment for the 
lower Mach numbers. For higher Mach numbers in the 
vicinity of the transonic dip region, the inviscid codes 
become increasingly conservative. For this wing, inviscid 
calculations should not be used for M > 0.80. 

2. For Mach numbers at and below the minimum 
transonic flutter speed index, the viscous methods, CAP- 
TSDV and CFL-3D are in agreement and both provide 
good agreement with experiment, largely correcting the 
deficiency in the inviscid methods. Also, the finite ele- 
ment structural model was not updated with information 
from the model vibration testing. This may account for a 
significant portion of the remaining differences between 
experiment and the calculations. 

3. Linear flutter calculations 43 are in excellent agree- 
ment with experiment up to M = 0.85, but cannot be 
relied upon for higher transonic Mach numbers. The 
good agreement in the lower transonic speed range is 
due to well-known compensating defects of linear theory 
wherein thickness and viscous effects are neglected. 
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Pressure Displacement Thickness Skin Friction 



Lower Surface 


Fig. 18 Contour plots of business jet wing pressure and boundary layer parameters at statically deformed conditions: M = 
0.888, Q = 79 psf n a « 0.2°, Rec = 1.14 million. 


All of the results discussed thus far were obtained 
from transient or harmonic responses of small ampli- 
tude, that is, wing tip response amplitudes were less 
than several tenths of an inch. Under these conditions, 
no large changes of the static aerodynamic loading oc- 
curred and transient responses exhibited exponential sta- 
bility, characteristic of a “locally linear” system behavior. 
At M = 0.888 the CAP-TSDV code was able to cal- 
culate large amplitude response motions which demon- 
strated limit cycle behavior. The motion was calculated 
for the experimental flutter dynamic pressure of 79 lb/ft 2 . 
The conditions for the limit cycle are noted in Fig. 19 
by the solid symbol indicating a 0.5 Hz. increase in fre- 
quency over the small amplitude value. Figure 20 shows 
two transient responses confirming the limit cycle behav- 
ior. The motions were excited from converged statically 
deformed conditions by multiplying the modal displace- 
ments and velocities by factors of 5.0 for Fig. 20a and 
0.5 for Fig. 20b. The larger factor simulates a wing 
tip displacement of about 7 inches, resulting in decay- 
ing oscillations to a limit cycle with an amplitude of 5—6 
inches peak-to-peak. The smaller factor results in os- 
cillations growing in amplitude to the limit cycle. This 
behavior is similar to model behavior observed during the 
test. Video tape of the model motions at the experimen- 
tal “flutter” conditions for this Mach number shows the 


model to be undergoing constant amplitude wing oscilla- 
tions with amplitude of slightly less than one tip chord 
(6.3 inches) peak-to-peak. This is in very good agree- 
ment with the calculated LCO amplitude and frequency 
shown in Fig. 20. The plate construction of the model 
provides sufficient strength to allow the model to sus- 
tain oscillations of this amplitude without structural fail- 
ure. Inspection of the wing boundary layer parameters 
and surface pressures during the calculated limit cycle 
oscillations confirmed that the flow over the wing was 
intermittently separating and reattaching in the outboard 
upper and lower surface regions described above. This 
apparently provides the mechanism needed to quench the 
growth of the unstable flutter mode motions. 
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(a) flutter speed index. 



(b) frequency. 

Fig. 19 Comparison between experimental and calculated 
flutter speed index and frequency for a business jet flutter 
model tested in air. 


Concluding Remarks 

A viscous-inviscid interactive coupling method has 
been described, directed towards the computation of un- 



time, sec. 

(a) amplitude decaying to limit cycle oscillation. 



time, sec. 


(b) amplitude growing to limit cycle oscillation. 

Fig. 20 Calculated limit cycle response for a business jet 
wing flutter model: M = 0.888, Q s 79 psf., Re* = 1.14 
million. 

steady separating and reattaching transonic flows which 
must be treated in cases of self-excited shock-induced os- 
cillations and transonic flutter. Lag-entrainment integral 
boundary layer equations and a transonic small distur- 
bance potential code are coupled with a variable gain, in- 
tegral control coupling method. The buffet onset bound- 
ary for the NACA 0012 airfoil is shown computationally 
to be a Hopf bifurcation and good agreement with exper- 
iment is shown for the boundary and oscillation frequen- 
cies. For the 18 percent thick circular arc airfoil, very 
good agreement with experiment is shown for shock os- 
cillation frequencies. The hysteresis with Mach number 
of the oscillation onset boundary is reproduced compu- 
tationally, including an amplitude threshold, jump phe- 
nomenon stability boundary. 

Rutter calculations for the AGARD 445.6 flutter 
model are in excellent agreement with experiment for 
M < 1.0 for models tested in air and heavy gas. Cal- 
culations with the CAP-TSDV code are in excellent 
agreement with results from a Navier-Stokes code at 
M = 0.96. For Mach numbers below and very near 
unity, viscous modeling is required for such thin wings 
in order to achieve acceptable accuracy. In this region, 
calculations show evidence of small amplitude limit cy- 
cle behavior. For very low supersonic Mach numbers, 
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agreement with experiment is not yet satisfactory. Im- 
proved modeling and/or knowledge of wind tunnel test 
conditions is needed. 

Flutter calculations for a business jet wing model 
also show very good agreement with experiment for the 
available test data up to M — 0.9. For this thicker 
wing, the requirement for viscous modeling extends to 
lower transonic Mach numbers. Again, calculations with 
the CAP-TSDV code are in very good agreement with 
a Navier-Stokes code at M = 0.888 for small amplitude 
flutter motions. For large amplitude wing oscillations, the 
CAP-TSDV code predicts limit cycle behavior in very 
good agreement with that observed during wind tunnel 
tests of the model. 
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